Introduction to the hands-on project:

Cognitive skills impact on fMRI restings state intrinsic connectivity

The aim of this 6-day project is to familiarize you with the basic steps of functional connectivity analysis by working with openly available tools and data.

Throughout the course example datasets will be used in these notebooks. However, all the notebooks can run the necessary code, so feel free to alter the dataset, analytic approach, visualization, etc... as much as you would like.

Project Schedule:

Day Topic Notebook
1 Downloading data 1_selecting_data.ipynb
1 Beginning preprocessing 2_preprocessing_mri_data.ipynb
2 Connectivity analysis 4_computing_connectivity.ipynb
3 Working with behavioral data 3_behavioral_data.ipynb
4-5 Group-level analyses 5_group-level_analysis.ipynb
6 Data visualization techniques 6_visualizing_results.ipynb

1. Selecting data

There is plenty of openly available MRI and behavioral data openly available for download.
To begin, using the Chrome browser, go to: openneuro.org

One of the advantages of all data available at openneuro.org is that it is already prepared in BIDS format, which is a standard directory structure and file naming convention that enables the easy use of a broad range of analysis tools called BIDS Apps. (We'll be using the BIDS App fmriprep for MRI data preprocessing in the next section.)

Choosing a dataset

When searching through openneuro.org, look for a dataset that has a behavioral measure that is of interest to you.

Available datasets at openneuro.org

Other resources for finding a BIDS-compatible dataset:

Downloading

Download the dataset, unzip it, and place it into the data directory.

In [ ]:
! datalad install -g https://github.com/OpenNeuroDatasets/ds000245.git

Preprocessing using fmriprep

fmriprep allows us to use optimized pipelines...

Open a terminal window and execute the following command:

In [ ]:
# Faster version using '--sloppy' (this is only used in the course due to our tight schedule): 
docker run -ti --rm \
  -v /home/padawan/cajalcourse/ds000245:/data:ro \
  -v /home/padawan/cajalcourse/ds000245_preproc:/out \
  -v /home/padawan/cajalcourse:/fs poldracklab/fmriprep:latest \
  --fs-no-reconall  \
  --output-spaces MNI152NLin2009cAsym \
  --force-no-bbr \
  --dummy-scans 4 \
  --sloppy \
  --notrack \
  --n_cpus 6 \
  --mem_mb 12000 \
  --participant_label sub-CTL01 sub-CTL02 sub-CTL03 sub-CTL04 \
  --fs-license-file /fs/license.txt \
  /data /out participant
In [ ]:
 
In [2]:
{
    "tags": [
        "hide_input",
    ]
}

a = 1
b = 3

print(a)
1
In [ ]:
 
In [ ]:
 

Quality control

fmriprep outputs various quality control metrics on the participant-level that are viewable through an html page.

Working with behavioral/phenotypic data

1. Loading the behavioral data

Locate the file containing the behavioral data .tsv or .csv file:

In [1]:
import pandas as pd


subject = 'CON02'
confounds_file = '/data/rw_eleves/CajalProject12/preproc/fmriprep/sub-%s/ses-postop/func/sub-%s_ses-postop_task-rest_desc-confounds_regressors.tsv' % (subject, subject)              
In [2]:
print(confounds_file)
/data/rw_eleves/CajalProject12/preproc/fmriprep/sub-CON02/ses-postop/func/sub-CON02_ses-postop_task-rest_desc-confounds_regressors.tsv
In [3]:
cf_all = pd.read_csv(confounds_file, delimiter='\t')
In [7]:
cf_all.columns
Out[7]:
Index(['csf', 'csf_derivative1', 'csf_power2', 'csf_derivative1_power2',
       'white_matter', 'white_matter_derivative1',
       'white_matter_derivative1_power2', 'white_matter_power2',
       'global_signal', 'global_signal_derivative1',
       ...
       'motion_outlier05', 'motion_outlier06', 'motion_outlier07',
       'motion_outlier08', 'motion_outlier09', 'motion_outlier10',
       'motion_outlier11', 'motion_outlier12', 'motion_outlier13',
       'motion_outlier14'],
      dtype='object', length=205)
In [5]:
cf = cf_all.filter(['csf',
                    'white_matter',
                    'global_signal', 
                    'trans_x', 
                    'trans_y', 
                    'trans_z', 
                    'rot_x', 
                    'rot_y', 
                    'rot_z'], axis=1)
In [8]:
cf.to_csv(path_or_buf='/data/rw_eleves/CajalProject12/preproc/fmriprep/sub-%s/ses-postop/func/sub-%s_ses-postop_task-rest_desc-confounds_regressors_selected.csv' % (subject,subject), index=False)
In [11]:
confounds_name = '/data/rw_eleves/CajalProject12/preproc/fmriprep/sub-%s/ses-postop/func/sub-%s_ses-postop_task-rest_desc-confounds_regressors_selected.csv' % (subject,subject)
In [12]:
confounds_name
Out[12]:
'/data/rw_eleves/CajalProject12/preproc/fmriprep/sub-CON02/ses-postop/func/sub-CON02_ses-postop_task-rest_desc-confounds_regressors_selected.csv'
In [13]:
import sklearn
import sklearn.datasets

# use load command from nilearn
# not glob, because it reorders

dataset = sklearn.datasets.base.Bunch(func=['/data/rw_eleves/CajalProject12/preproc/fmriprep/sub-%s/ses-postop/func/sub-%s_ses-postop_task-rest_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz' % (subject, subject)], 
                                      confounds=[confounds_name])
In [94]:
import sklearn
import sklearn.datasets

# use load command from nilearn
# not glob, because it reorders

dataset = sklearn.datasets.base.Bunch(func=['/data/rw_eleves/CajalProject12/ABIDE_pcp/cpac/nofilt_noglobal/Olin_0050102_func_preproc.nii.gz'], 
                                      confounds=['BABAB'])
In [95]:
dataset
Out[95]:
{'func': ['/data/rw_eleves/CajalProject12/ABIDE_pcp/cpac/nofilt_noglobal/Olin_0050102_func_preproc.nii.gz'],
 'confounds': ['BABAB']}
In [42]:
%matplotlib inline

import nilearn
from nilearn import plotting
import joblib
from nilearn import input_data
import numpy as np

motor_coords = [(-36, -28, 56)]
pcc_coords = [(0, -52, 18)]
ba45_coords = [(-52, 28, 2)]
precun = [(0, -66, 46)]

seed_masker = input_data.NiftiSpheresMasker(
    pcc_coords, radius=4,
    detrend=True, standardize=True,
    low_pass=0.1, high_pass=0.01, t_r=2.4,
    memory='nilearn_cache', memory_level=1, verbose=0)

seed_time_series = seed_masker.fit_transform(dataset.func[0], confounds=[dataset.confounds[0]])

brain_masker = input_data.NiftiMasker(
    smoothing_fwhm=6,
    detrend=True, standardize=True,
    low_pass=0.1, high_pass=0.01, t_r=2.4,
    memory='nilearn_cache', memory_level=1, verbose=0)

brain_time_series = brain_masker.fit_transform(dataset.func[0], confounds=[dataset.confounds[0]]) 


seed_to_voxel_correlations = (np.dot(brain_time_series.T, seed_time_series) /
                              seed_time_series.shape[0]
                              )

seed_to_voxel_correlations_img = brain_masker.inverse_transform(
    seed_to_voxel_correlations.T)
display = plotting.plot_stat_map(seed_to_voxel_correlations_img,
                                 threshold=0.5, vmax=1,
                                 cut_coords=pcc_coords[0],
                                 title="Seed-to-voxel correlation (pcc seed)"
                                 )

display.add_markers(marker_coords=pcc_coords, marker_color='g',
                    marker_size=70)

# At last, we save the plot as pdf.
display.savefig('pcc_coords_seed_correlation.pdf')
In [107]:
#Copy of the before steps for the ABIDE data
import nilearn
from nilearn import plotting
import joblib
from nilearn import input_data
import numpy as np

motor_coords = [(-36, -28, 56)]
pcc_coords = [(0, -52, 18)]
ba45_coords = [(-52, 28, 2)]
precun = [(0, -66, 46)]
L_DMN = [(-45.8, -64.78,31.84)]

seed_masker = input_data.NiftiSpheresMasker(
   L_DMN, radius=4,
    detrend=True, standardize=True,
    low_pass=0.1, high_pass=0.01, t_r=2,
    memory='nilearn_cache', memory_level=1, verbose=0)

seed_time_series = seed_masker.fit_transform(dataset.func[0] )

brain_masker = input_data.NiftiMasker(
    smoothing_fwhm=6,
    detrend=True, standardize=True,
    low_pass=0.1, high_pass=0.01, t_r=2,
    memory='nilearn_cache', memory_level=1, verbose=0)

brain_time_series = brain_masker.fit_transform(dataset.func[0]) 


seed_to_voxel_correlations = (np.dot(brain_time_series.T, seed_time_series) /
                              seed_time_series.shape[0]
                              )

seed_to_voxel_correlations_img = brain_masker.inverse_transform(
    seed_to_voxel_correlations.T)
display = plotting.plot_stat_map(seed_to_voxel_correlations_img,
                                 threshold=0.4, vmax=1,
                                 cut_coords= L_DMN[0],
                                 title="Seed-to-voxel correlation (L DMN seed)"
                                 )

display.add_markers(marker_coords=pcc_coords, marker_color='g',
                    marker_size=70)

# At last, we save the plot as pdf.
display.savefig('L_DMN_coords_seed_correlation_ABIDE_Test.pdf')
In [104]:
display = plotting.plot_stat_map(seed_to_voxel_correlations_img,
                                 threshold=0.8, vmax=1,
                                 cut_coords=pcc_coords[0],
                                 title="Seed-to-voxel correlation (pcc seed)"
                                 )

display.add_markers(marker_coords=pcc_coords, marker_color='g',
                    marker_size=70)

# At last, we save the plot as pdf.
display.savefig('pcc_coords_seed_correlation_ABIDE_Test.pdf')
In [73]:
# function to define a seed using its MNI coordinates and return voxel-wise correlations 
# with the seed, and save the image

def seed_cor(coords,name):
    seed_masker = input_data.NiftiSpheresMasker(
    coords, radius=4,
    detrend=True, standardize=True,
    low_pass=0.1, high_pass=0.01, t_r=2.4,
    memory='nilearn_cache', memory_level=1, verbose=0)

    seed_time_series = seed_masker.fit_transform(dataset.func[0], confounds=[dataset.confounds[0]])

    brain_masker = input_data.NiftiMasker(
    smoothing_fwhm=6,
    detrend=True, standardize=True,
    low_pass=0.1, high_pass=0.01, t_r=2.4,
    memory='nilearn_cache', memory_level=1, verbose=0)

    brain_time_series = brain_masker.fit_transform(dataset.func[0], confounds=[dataset.confounds[0]]) 


    seed_to_voxel_correlations = (np.dot(brain_time_series.T, seed_time_series) /
                              seed_time_series.shape[0]
                              )

    seed_to_voxel_correlations_img = brain_masker.inverse_transform(
    seed_to_voxel_correlations.T)
    display = plotting.plot_stat_map(seed_to_voxel_correlations_img,
                                 threshold=0.5, vmax=1,
                                 cut_coords=coords[0],
                                 title="Seed-to-voxel correlation %s" % name
                                 )

    display.add_markers(marker_coords=coords, marker_color='g',
                    marker_size=70)

# At last, we save the plot as pdf.
    display.savefig('%s_coords_seed_correlation.pdf' % name)

    return display
In [84]:
seed_cor(coords=[(-20, -20, 20)],name = 'random')
Out[84]:
<nilearn.plotting.displays.OrthoSlicer at 0x7eff90addac8>
In [111]:
seed_to_voxel_correlations_fisher_z = np.arctanh(seed_to_voxel_correlations)
print("Seed-to-voxel correlation Fisher-z transformed: min = %.3f; max = %.3f"
      % (seed_to_voxel_correlations_fisher_z.min(),
         seed_to_voxel_correlations_fisher_z.max()
         )
      )

# Finally, we can tranform the correlation array back to a Nifti image
# object, that we can save.
seed_to_voxel_correlations_fisher_z_img = brain_masker.inverse_transform(
    seed_to_voxel_correlations_fisher_z.T)
seed_to_voxel_correlations_fisher_z_img.to_filename(
    'L_DMN_ABIDE_seed_correlation_z.nii.gz') 

view = plotting.view_img_on_surf(seed_to_voxel_correlations_fisher_z_img, 
                                 threshold='80%')
view
Seed-to-voxel correlation Fisher-z transformed: min = -0.576; max = 1.641
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/nilearn/plotting/js_plotting_utils.py:103: UserWarning: It seems you have created more than 10 nilearn views. As each view uses dozens of megabytes of RAM, you might want to delete some of them.
  warnings.warn('It seems you have created more than 10 '
Out[111]:

Day 3

Starting to develop our connectivity matrices for the ABIDE data set

In [4]:
from nilearn import datasets
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you may need to re-serialize those models with scikit-learn 0.21+.
  warnings.warn(msg, category=DeprecationWarning)

Compute connectivity matrix on single dataset

In [222]:
# Fetch all subjects 

ABIDE_data = datasets.fetch_abide_pcp(data_dir='/data/rw_eleves/CajalProject12', pipeline='cpac', 
                                 n_subjects = None, band_pass_filtering=False, global_signal_regression=False, derivatives=['func_preproc'], 
                                 quality_checked=True, url=None, verbose=0)
In [223]:
# Fetch the ROI-atlas (Harvard Oxford)

atlas = datasets.fetch_atlas_msdl()
atlas_filename = atlas.maps
labels = atlas['labels']
In [228]:
# Defining ROIs as seeds
from nilearn.input_data import NiftiMapsMasker

masker = NiftiMapsMasker(maps_img=atlas_filename, standardize=True, smoothing_fwhm=6,
    detrend=True, low_pass=0.1, high_pass=0.01, t_r=2,
    memory='nilearn_cache', verbose=5)

time_series = masker.fit_transform(ABIDE_data.func_preproc[0])
[NiftiMapsMasker.fit_transform] loading regions from /home/padawan/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii
Resampling maps
[Memory]0.0s, 0.0min    : Loading resample_img...
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/nilearn/_utils/cache_mixin.py:302: UserWarning: memory_level is currently set to 0 but a Memory object has been provided. Setting memory_level to 1.
  warnings.warn("memory_level is currently set to 0 but "
________________________________________resample_img cache loaded - 0.6s, 0.0min
[Memory]0.8s, 0.0min    : Loading filter_and_extract...
__________________________________filter_and_extract cache loaded - 0.0s, 0.0min

Only for n=1

In [35]:
# Calculate correlation example for n=1
from nilearn.connectome import ConnectivityMeasure
correlation_measure = ConnectivityMeasure(kind='correlation')
correlation_matrix = correlation_measure.fit_transform([time_series])[0]

# Define positive matrix
correlation_matrix_pos = correlation_matrix.copy()
correlation_matrix_pos[correlation_matrix < 0.0] = 0.0

# Define negative matrix
correlation_matrix_neg = correlation_matrix.copy()
correlation_matrix_neg[correlation_matrix > 0.0] = 0.0

# Display the full correlation matrix
import numpy as np
from nilearn import plotting
# Mask out the major diagonal
np.fill_diagonal(correlation_matrix, 0)
plotting.plot_matrix(correlation_matrix, labels=labels, colorbar=True,
                     vmax=0.8, vmin=-0.8)

# Display the positive correlation matrix
np.fill_diagonal(correlation_matrix_pos, 0)
plotting.plot_matrix(correlation_matrix_pos, labels=labels, colorbar=True,
                     vmax=0.8, vmin=-0.8)


# Display the negative correlation matrix
np.fill_diagonal(correlation_matrix_neg, 0)
plotting.plot_matrix(correlation_matrix_neg, labels=labels, colorbar=True,
                     vmax=0.8, vmin=-0.8)
Out[35]:
<matplotlib.image.AxesImage at 0x7f8feb3b35f8>
In [1]:
# Plotting the connectome
from nilearn import plotting
coords = atlas.region_coords

# We threshold to keep only the 20% of edges with the highest value
# because the graph is very dense
plotting.plot_connectome(correlation_matrix_pos, coords,
                         edge_threshold="90%", colorbar=True)

plotting.show()

# We threshold to keep only the 20% of edges with the highest value
# because the graph is very dense
plotting.plot_connectome(correlation_matrix_neg, coords,
                         edge_threshold="90%", colorbar=True)

plotting.show()
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you may need to re-serialize those models with scikit-learn 0.21+.
  warnings.warn(msg, category=DeprecationWarning)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-087704398e7e> in <module>
      7 # Plotting the connectome
      8 from nilearn import plotting
----> 9 coords = atlas.region_coords
     10 
     11 # We threshold to keep only the 20% of edges with the highest value

NameError: name 'atlas' is not defined

3. Computing functional connectivity matrices

Now that we have our preprocessed data, we can begin the functional connectivity analysis.

We begin by loading the data:

In [229]:
# Matrix plotting from Nilearn: nilearn.plotting.plot_matrix
import numpy as np
import matplotlib.pylab as plt

def plot_matrices(matrices, matrix_kind):
    n_matrices = len(matrices)
    fig = plt.figure(figsize=(n_matrices * 4, 4))
    for n_subject, matrix in enumerate(matrices):
        plt.subplot(1, n_matrices, n_subject + 1)
        matrix = matrix.copy()  # avoid side effects
        # Set diagonal to zero, for better visualization
        np.fill_diagonal(matrix, 0)
        vmax = np.max(np.abs(matrix))
        title = '{0}, subject {1}'.format(matrix_kind, n_subject)
        plotting.plot_matrix(matrix, vmin=-vmax, vmax=vmax, cmap='RdBu_r',
                             title=title, figure=fig, colorbar=True)
In [230]:
msdl_data = datasets.fetch_atlas_msdl()
msdl_coords = msdl_data.region_coords
n_regions = len(msdl_coords)
print('MSDL has {0} ROIs, part of the following networks :\n{1}.'.format(
    n_regions, msdl_data.networks))
MSDL has 39 ROIs, part of the following networks :
[b'Aud', b'Aud', b'Striate', b'DMN', b'DMN', b'DMN', b'DMN', b'Occ post', b'Motor', b'R V Att', b'R V Att', b'R V Att', b'R V Att', b'Basal', b'L V Att', b'L V Att', b'L V Att', b'D Att', b'D Att', b'Vis Sec', b'Vis Sec', b'Vis Sec', b'Salience', b'Salience', b'Salience', b'Temporal', b'Temporal', b'Language', b'Language', b'Language', b'Language', b'Language', b'Cereb', b'Dors PCC', b'Cing-Ins', b'Cing-Ins', b'Cing-Ins', b'Ant IPS', b'Ant IPS'].
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/numpy/lib/npyio.py:2322: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.
  output = genfromtxt(fname, **kwargs)
In [231]:
from nilearn import input_data

masker = input_data.NiftiMapsMasker(
    msdl_data.maps, resampling_target="data", t_r=2, detrend=True,
    low_pass=.1, high_pass=.01, memory='nilearn_cache', memory_level=1)
In [232]:
# Pulling out the groups + nec info
ABIDE_subjects = []
pooled_subjects = []
control_subjects = []
ABIDE_labels = []  # 1 if Autism, 2 if control
for func_file,  phenotypic in zip(
        ABIDE_data.func_preproc, ABIDE_data.phenotypic):
    time_series = masker.fit_transform(func_file)
    pooled_subjects.append(time_series)
    is_ABIDE = phenotypic['DX_GROUP']
    if is_ABIDE == 1:
        ABIDE_subjects.append(time_series)
    else:
        control_subjects.append(time_series)
    
    ABIDE_labels.append(is_ABIDE)

print('Data has {0} ABIDE subjects.'.format(len(ABIDE_subjects)))
print('Data has {0} control subjects.'.format(len(control_subjects)))
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/nilearn/input_data/nifti_maps_masker.py:326: UserWarning: Persisting input arguments took 0.59s to run.
If this happens often in your code, it can cause performance problems 
(results will be correct in all cases). 
The reason for this is probably some large input arguments for a wrapped
 function (e.g. large strings).
THIS IS A JOBLIB ISSUE. If you can, kindly provide the joblib's team with an
 example so that they can fix the problem.
  verbose=self.verbose)
Data has 320 ABIDE subjects.
Data has 428 control subjects.
In [307]:
#Had one extra file that we weren't using 
#this made a list for all the files we had

SUBS = []
for i in range(len(ABIDE_data.phenotypic)):
        SUBS.append(ABIDE_data.phenotypic[i][2])
        
#Compares the ID numbers in SUBS with the ID numbers in Plist to identify the extra file

print(set(Plist).difference(set(SUBS)))
{50026}
In [233]:
from nilearn.connectome import ConnectivityMeasure
In [236]:
#creating correlation matrices for Autism subjects

correlation_measure_ABIDE = ConnectivityMeasure(kind='correlation')
correlation_matrices_ABIDE = correlation_measure_ABIDE.fit_transform(ABIDE_subjects)

# All individual coefficients are stacked in a unique 2D matrix.
print('Correlations of ABIDE patients are stacked in an array of shape {0}'
      .format(correlation_matrices_ABIDE.shape))

#Take mean correlation across autism group

mean_correlation_matrix_ABIDE = correlation_measure_ABIDE.mean_
print('Mean correlation has shape {0}.'.format(mean_correlation_matrix_ABIDE.shape))

from nilearn import plotting

#Plot the first 4 autism subjects and the the mean connectome across the autism group
plot_matrices(correlation_matrices_ABIDE[:4], 'correlation')
plotting.plot_connectome(mean_correlation_matrix_ABIDE, msdl_coords, edge_threshold = 0.6,
                         title='mean correlation over 320 ABIDE subjects')

mean_correlation_matrix_extra_ABIDE = np.expand_dims(mean_correlation_matrix_ABIDE, axis=0)
plot_matrices(mean_correlation_matrix_extra_ABIDE, 'correlation')
Correlations of ABIDE patients are stacked in an array of shape (320, 39, 39)
Mean correlation has shape (39, 39).
In [237]:
#creating correlation matrices for control subjects
correlation_measure_control = ConnectivityMeasure(kind='correlation')
correlation_matrices_control = correlation_measure_control.fit_transform(control_subjects)

# All individual coefficients are stacked in a unique 2D matrix.
print('Correlations of control patients are stacked in an array of shape {0}'
      .format(correlation_matrices_control.shape))

#Take mean correlation across control group

mean_correlation_matrix_control = correlation_measure_control.mean_
print('Mean correlation has shape {0}.'.format(mean_correlation_matrix_control.shape))

#Plot the first 4 control subjects and the the mean connectome across the control group
plot_matrices(correlation_matrices_control[:4], 'correlation')
plotting.plot_connectome(mean_correlation_matrix_control, msdl_coords, edge_threshold = 0.6,
                         title='mean correlation over 428 control subjects')

mean_correlation_matrix_extra_control = np.expand_dims(mean_correlation_matrix_control, axis=0)
plot_matrices(mean_correlation_matrix_extra_control, 'correlation')
Correlations of control patients are stacked in an array of shape (428, 39, 39)
Mean correlation has shape (39, 39).

---Day 4---

Partial and Tangent correlation for both groups + plots

In [240]:
#partial and tangent correlation for the Autism group

partial_correlation_measure_ABIDE = ConnectivityMeasure(kind='partial correlation')

partial_correlation_matrices_ABIDE = partial_correlation_measure_ABIDE.fit_transform(
    ABIDE_subjects)

plot_matrices(partial_correlation_matrices_ABIDE[:4], 'partial')
plotting.plot_connectome(
    partial_correlation_measure_ABIDE.mean_, msdl_coords,
    title='mean partial correlation over 320 ABIDE subjects')

tangent_measure_ABIDE = ConnectivityMeasure(kind='tangent')
tangent_matrices_ABIDE = tangent_measure_ABIDE.fit_transform(ABIDE_subjects)
plot_matrices(tangent_matrices_ABIDE[:4], 'tangent variability')
plotting.plot_connectome(
    tangent_measure_ABIDE.mean_, msdl_coords,
    title='mean tangent connectivity over 320 ABIDE subjects')
Out[240]:
<nilearn.plotting.displays.OrthoProjector at 0x7f8fcaa6e978>
In [241]:
#partial and tangent correlation for the Control group

partial_correlation_measure_control = ConnectivityMeasure(kind='partial correlation')

partial_correlation_matrices_control = partial_correlation_measure_control.fit_transform(
    control_subjects)

plot_matrices(partial_correlation_matrices_control[:4], 'partial')
plotting.plot_connectome(
    partial_correlation_measure_control.mean_, msdl_coords,
    title='mean partial correlation over 428 control subjects')

tangent_measure_control = ConnectivityMeasure(kind='tangent')
tangent_matrices_control = tangent_measure_control.fit_transform(control_subjects)
plot_matrices(tangent_matrices_control[:4], 'tangent variability')
plotting.plot_connectome(
    tangent_measure_control.mean_, msdl_coords,
    title='mean tangent connectivity over 428 control subjects')
Out[241]:
<nilearn.plotting.displays.OrthoProjector at 0x7f8fb7766278>
In [242]:
connectivity_biomarkers = {}
kinds = ['correlation', 'partial correlation', 'tangent']
for kind in kinds:
    conn_measure = ConnectivityMeasure(kind=kind, vectorize=True)
    connectivity_biomarkers[kind] = conn_measure.fit_transform(pooled_subjects)

# For each kind, all individual coefficients are stacked in a unique 2D matrix.
print('{0} correlation biomarkers for each subject.'.format(
    connectivity_biomarkers['correlation'].shape[1]))
780 correlation biomarkers for each subject.
In [243]:
from sklearn.model_selection import StratifiedKFold

"""
classes = ['{0}{1}'.format(ABIDE_label)
           for ABIDE_label in zip(ABIDE_labels)]
"""

cv = StratifiedKFold(n_splits=3)
In [244]:
# Note that in cv.split(X, y),
# providing y is sufficient to generate the splits and
# hence np.zeros(n_samples) may be used as a placeholder for X
# instead of actual training data.
from sklearn.svm import LinearSVC
from sklearn.model_selection import cross_val_score

mean_scores = []
for kind in kinds:
    svc = LinearSVC(random_state=0)
    cv_scores = cross_val_score(svc,
                                connectivity_biomarkers[kind],
                                y=ABIDE_labels,
                                cv=cv,
                                groups=ABIDE_labels,
                                scoring='accuracy',
                                )
    mean_scores.append(cv_scores.mean())
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/sklearn/svm/base.py:929: ConvergenceWarning: Liblinear failed to converge, increase the number of iterations.
  "the number of iterations.", ConvergenceWarning)
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/sklearn/svm/base.py:929: ConvergenceWarning: Liblinear failed to converge, increase the number of iterations.
  "the number of iterations.", ConvergenceWarning)
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/sklearn/svm/base.py:929: ConvergenceWarning: Liblinear failed to converge, increase the number of iterations.
  "the number of iterations.", ConvergenceWarning)
In [245]:
from nilearn.plotting import show

plt.figure(figsize=(6, 4))
positions = np.arange(len(kinds)) * .1 + .1
plt.barh(positions, mean_scores, align='center', height=.05)
yticks = [kind.replace(' ', '\n') for kind in kinds]
plt.yticks(positions, yticks)
plt.xlabel('Classification accuracy')
plt.grid(True)
plt.tight_layout()

show()
In [253]:
# We use glob and split to only get the Sub_IDs in the data that we have (i.e. those that passed QC)

import glob,os

# Use glob to get the names from the folder
path = '/data/rw_eleves/CajalProject12/ABIDE_pcp/cpac/nofilt_noglobal'

files = [f for f in glob.glob(path + "**/*.nii.gz", recursive=True)]

# Get and sort subject ID
Plist=[]
for test in files:
    base=os.path.basename(test)
    basenames = base.split('_')
    if (basenames[1]=="a") or (basenames[1]=="b") or (basenames[1]=="c") or (basenames[1]=="d") or (basenames[1]=="1") or (basenames[1]=="2"):
        Plist.append(int(basenames[2]))
    else:
        Plist.append(int(basenames[1]))

Plist = list(np.sort(Plist))

Working with behavioral/phenotypic data

1. Loading the behavioral data

In [308]:
# Locate the file containing the behavioral data `.tsv` or `.csv` file:

import pandas as pd

behav = pd.read_csv('/data/rw_eleves/CajalProject12/ABIDE_pcp/Phenotypic_V1_0b_preprocessed1.csv', delimiter=',', na_values="-9999")

behav_short = behav.filter(['SUB_ID',
                   'subject', 'SITE_ID','DX_GROUP', 'DSM_IV_TR', 'AGE_AT_SCAN',  'SEX',
             'HANDEDNESS_CATEGORY', 'FIQ','VIQ', 'PIQ', 'FIQ_TEST_TYPE', 'VIQ_TEST_TYPE',
             'PIQ_TEST_TYPE',], axis=1)

behav_short = behav_short[behav_short['SUB_ID'].isin(Plist)] 

2. Check the distributions

In [309]:
# Locate the measure-of-interest:

%matplotlib inline

# Next, plot a histogram of the values:

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

f, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(18,5))
sns.violinplot(x=behav_short.DX_GROUP,y=behav_short.FIQ, ax=ax1)
sns.violinplot(x=behav_short.DX_GROUP,y=behav_short.VIQ,ax=ax2)
sns.violinplot(x=behav_short.DX_GROUP,y=behav_short.PIQ,ax=ax3)

fg, (ax4, ax5, ax6) = plt.subplots(1,3, figsize=(18,5))
sns.catplot(x="DX_GROUP",y="FIQ",data=behav_short, ax=ax4)
sns.catplot(x="DX_GROUP",y="VIQ",data=behav_short,ax=ax5)
sns.catplot(x="DX_GROUP",y="PIQ",data=behav_short,ax=ax6)
plt.close(4)
plt.close(5)
plt.close(3)

# needs to be fixed (to go on top)
'''
fh, (ax4, ax5, ax6) = plt.subplots(1,3, figsize=(18,5))
sns.distplot(behav_short.FIQ)
plt.show()
sns.distplot(behav_short.VIQ)
plt.show()
sns.distplot(behav_short.PIQ)
plt.show()
'''

f.tight_layout()
fg.tight_layout()


K, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(18,5))

sns.countplot(x='SEX', data=behav_short, ax=ax1)
sns.distplot(behav.AGE_AT_SCAN, ax=ax2)
sns.countplot(x='HANDEDNESS_CATEGORY',data=behav_short, ax=ax3 )

K.tight_layout()
plt.show()

Actual RS + behav

Q1: Relationship between IQ and

Does IQ influence resting state connectivity in people with/without ASD differently?

Similarity matrix --> modularity index

In [311]:
from community import community_louvain
import networkx as nx
In [312]:
part = community_louvain.best_partition(G)


mod = community_louvain.modularity(part,G)
print("modularity:", mod)
modularity: 0.8055540757328388
In [313]:
def modfunction(data): 
    Modscores=[] 

    for person in data:
        G = nx.from_numpy_matrix(person)
        part = community_louvain.best_partition(G)
        mod = community_louvain.modularity(part,G)
        Modscores.append(mod)    
    return Modscores
In [314]:
Modscore_Autism = modfunction(partial_correlation_matrices_ABIDE)

Modscore_control = modfunction(partial_correlation_matrices_control)

sns.distplot(Modscore_Autism)

sns.distplot(Modscore_control)
Out[314]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f8fbef40b38>
In [319]:
from scipy import stats

# Statistics group-level
stats.ttest_ind(Modscore_Autism,Modscore_control)
Out[319]:
Ttest_indResult(statistic=-1.7337800326727013, pvalue=0.08337026533289069)

Single Value Decompensation

In [315]:
import sklearn.decomposition
In [316]:
pca_outputs = []

for person in pooled_subjects:
    pca = sklearn.decomposition.PCA()
    pca.fit(person)
    pca_outputs.append(pca.explained_variance_ratio_[1])
        
In [320]:
Autism_pca = [pca_outputs[i] for i in np.where(behav_short.DX_GROUP==1)[0]]
Control_pca =  [pca_outputs[i] for i in np.where(behav_short.DX_GROUP==2)[0]]
stats.ttest_ind(Autism_pca,Control_pca)
Out[320]:
Ttest_indResult(statistic=-0.8375940124061025, pvalue=0.402527015111141)
In [321]:
sns.distplot(Autism_pca)
sns.distplot(Control_pca)
Out[321]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f8fb57f3e10>

Amplitude of Low Frequency Fluctuations

In [322]:
from scipy import signal
In [333]:
from scipy.integrate import simps
low, high = 0, 0.05
sf = 0.5 
win = 100 * sf

# we should refetch the data without tha bandpass filter
pooled_totalpower = []
for person in pooled_subjects:
    per_region = []
    for region in person.T:           
        freqs, psd = signal.welch(region, sf, nperseg=win)
        freq_res = freqs[5] - freqs[0]
        idx_testband = np.logical_and(freqs >= low, freqs <= high)
        testband_power = simps(psd[idx_testband], dx=freq_res)
        per_region.append(testband_power)
        
    pooled_totalpower.append(per_region)
In [342]:
totalpower = np.array(pooled_totalpower)
In [362]:
region_ttest = np.zeros((39,2))
    
for counter, region in enumerate(totalpower.T):
    Autism_power = [region[i] for i in np.where(behav_short.DX_GROUP==1)[0]]
    Control_power =  [region[i] for i in np.where(behav_short.DX_GROUP==2)[0]]
    t,p = stats.ttest_ind(Autism_power,Control_power)
            
    region_ttest[counter, 0] = t
    region_ttest[counter, 1] = p
    
In [363]:
region_ttest
Out[363]:
array([[-9.74424431e-02,  9.22401214e-01],
       [-6.42131264e-02,  9.48817728e-01],
       [ 1.22515433e-01,  9.02523821e-01],
       [-1.76431600e+00,  7.80879981e-02],
       [-3.01177641e+00,  2.68497857e-03],
       [-1.89510958e+00,  5.84642509e-02],
       [ 1.02962175e+00,  3.03521305e-01],
       [ 5.58869242e-01,  5.76418680e-01],
       [-2.15230405e+00,  3.16935594e-02],
       [ 8.10731816e-01,  4.17778260e-01],
       [ 3.17208995e+00,  1.57545710e-03],
       [ 3.01262232e+00,  2.67759743e-03],
       [ 6.59938520e-01,  5.09496976e-01],
       [-5.12080495e-01,  6.08746235e-01],
       [ 1.86495309e+00,  6.25804109e-02],
       [ 8.69088527e-01,  3.85078233e-01],
       [ 3.07934700e+00,  2.15063623e-03],
       [-8.02763071e-01,  4.22367364e-01],
       [-1.13999047e+00,  2.54656148e-01],
       [-2.56709563e-01,  7.97473769e-01],
       [-4.56456045e-01,  6.48194935e-01],
       [-2.32961671e-01,  8.15855072e-01],
       [ 2.19397909e+00,  2.85446434e-02],
       [-9.04924367e-02,  9.27920195e-01],
       [ 1.88564472e+00,  5.97311641e-02],
       [-5.65055589e-01,  5.72205835e-01],
       [-6.69649122e-01,  5.03288699e-01],
       [ 5.68227150e-01,  5.70051730e-01],
       [ 1.14873038e+00,  2.51035361e-01],
       [-2.01055949e-01,  8.40709631e-01],
       [ 6.12189519e-01,  5.40598917e-01],
       [ 1.20406379e+00,  2.28946822e-01],
       [-4.15104759e-01,  6.78184571e-01],
       [ 5.02710004e-02,  9.59919889e-01],
       [-1.88954124e-01,  8.50180166e-01],
       [-1.96692859e+00,  4.95616827e-02],
       [-2.11680677e+00,  3.46067121e-02],
       [-1.67816118e-01,  8.66773425e-01],
       [-8.24800661e-01,  4.09748420e-01]])
In [371]:
sig_ps = np.where(region_ttest[:,1] < 0.05)[0]
In [375]:
from nilearn import plotting, image

# First plot the map for the PCC: index 4 in the atlas
display = plotting.plot_stat_map(image.index_img(atlas_filename, sig_ps[0]),
                                 colorbar=False,
                                 title="ALFF Results in MSDL atlas")

# Now add as an overlay for other regions showing sig differences

for region in sig_ps[1:]:
    display.add_overlay(image.index_img(atlas_filename, region))

plotting.show()
    
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/numpy/ma/core.py:2788: UserWarning: Warning: converting a masked element to nan.
  order=order, subok=True, ndmin=ndmin)
In [379]:
masker = input_data.NiftiMapsMasker(atlas_filename, resampling_target="data")
In [383]:
import nibabel as nib
atlas_vol = nib.load(atlas_filename)
In [390]:
data = np.argmax(atlas_vol.get_fdata(),axis=3)

img = nib.Nifti1Image(data, atlas_vol.affine)
img.to_filename("atlas_vol.nii.gz")
In [418]:
data_res = np.asarray(data.copy(),dtype = float)


for i in range(39):
    data_res[data == i] = region_ttest[i,0]
        
img = nib.Nifti1Image(data_res, atlas_vol.affine)   
plotting.plot_stat_map(img, threshold = 1)
Out[418]:
<nilearn.plotting.displays.OrthoSlicer at 0x7f8fc27c65c0>
In [425]:
fsaverage = datasets.fetch_surf_fsaverage()

from nilearn import surface

texture = surface.vol_to_surf(img, fsaverage.pial_right)

view = plotting.view_img_on_surf(img, threshold=1)

view
Out[425]:
In [ ]:
 
In [337]:
'''# Plot the power spectrum
sns.set(font_scale=1.2, style='white')

plt.figure(figsize=(8, 4))
plt.plot(freqs, psd, color='y', lw=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (V^2 / Hz)')
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
plt.xlim([0, freqs.max()])
sns.despine()'''
Out[337]:
'# Plot the power spectrum\nsns.set(font_scale=1.2, style=\'white\')\nplt.figure(figsize=(8, 4))\nplt.plot(freqs, psd, color=\'y\', lw=2)\nplt.xlabel(\'Frequency (Hz)\')\nplt.ylabel(\'Power spectral density (V^2 / Hz)\')\nplt.ylim([0, psd.max() * 1.1])\nplt.title("Welch\'s periodogram")\nplt.xlim([0, freqs.max()])\nsns.despine()'
In [145]:
#Try again on one of regional time series


masker = NiftiMapsMasker(maps_img=atlas_filename, standardize=True, 
                         detrend=True, t_r=2, memory='nilearn_cache', verbose=5)

time_series_by_region = masker.fit_transform(ABIDE_data_for_ALLF.func_preproc[0])
[NiftiMapsMasker.fit_transform] loading regions from /home/padawan/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii
Resampling maps
[Memory]0.0s, 0.0min    : Loading resample_img...
/home/padawan/miniconda3/envs/cajal/lib/python3.7/site-packages/nilearn/_utils/cache_mixin.py:302: UserWarning: memory_level is currently set to 0 but a Memory object has been provided. Setting memory_level to 1.
  warnings.warn("memory_level is currently set to 0 but "
________________________________________resample_img cache loaded - 0.8s, 0.0min
________________________________________________________________________________
[Memory] Calling nilearn.input_data.base_masker.filter_and_extract...
filter_and_extract('/data/rw_eleves/CajalProject12/ABIDE_pcp/cpac/nofilt_noglobal/Pitt_0050003_func_preproc.nii.gz', 
<nilearn.input_data.nifti_maps_masker._ExtractionFunctor object at 0x7f974f76ff28>, 
{ 'allow_overlap': True,
  'detrend': True,
  'dtype': None,
  'high_pass': None,
  'low_pass': None,
  'maps_img': '/home/padawan/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii',
  'mask_img': None,
  'smoothing_fwhm': None,
  'standardize': True,
  't_r': 2,
  'target_affine': None,
  'target_shape': None}, confounds=None, dtype=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=5)
[NiftiMapsMasker.transform_single_imgs] Loading data from /data/rw_eleves/CajalProject12/ABIDE_pcp/cpac/nofilt_noglobal/Pitt_0050003_func_preproc.nii.gz
[NiftiMapsMasker.transform_single_imgs] Extracting region signals
[NiftiMapsMasker.transform_single_imgs] Cleaning extracted signals
_______________________________________________filter_and_extract - 7.8s, 0.1min
In [149]:
sf = 0.5  
time = np.shape(brain_time_seriesfor_ALFF)[0]/sf

win = 100 * sf
freqs, psd = signal.welch(time_series_by_region[:,30], sf, nperseg=win)

# Plot the power spectrum
sns.set(font_scale=1.2, style='white')
plt.figure(figsize=(8, 4))
plt.plot(freqs, psd, color='y', lw=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (V^2 / Hz)')
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
plt.xlim([0, freqs.max()])
sns.despine()
In [154]:
from scipy.integrate import simps

# Define delta lower and upper limits of the frequencies
low, high = 0, 0.05

# Find intersecting values in frequency vector
idx_testband = np.logical_and(freqs >= low, freqs <= high)

freq_res = freqs[5] - freqs[0] 


# Compute the absolute power by approximating the area under the curve
testband_power = simps(psd[idx_testband], dx=freq_res)
print('Absolute testband power: %.3f uV^2' % testband_power)
Absolute testband power: 1.283 uV^2